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We show how to describe the coupling of electrons to non-uniform magnetic fields in the frame- 
work of the widely used norm-conserving pseudopotential approximation for electronic structure 
calculations. Our derivation applies to magnetic fields that are smooth on the scale of the core 
region. The method is validated by application to the calculation of the magnetic susceptibility of 
molecules. Our results are compared with high quality all electron quantum chemical results, and 
another recently proposed formalism. 
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The coupling of the electrons in matter with prob- 
ing electromagnetic fields or charged particles provides 
the basis for nearly all known analytical experimental 
techniques. However, the most computationally efficient 
schemes for the first principles prediction (and hence in- 
terpretation) of the experimental observables require the 
use of approximate, and crucially, nonlocal Hamiltonians. 
Most notably, the use of nonlocal pseudopotentials, along 
with with density functional theory, is often referred to 
as the Standard Model of modern electronic structure 
theory. This is not without justification. The computa- 
tional efhciency, and accuracy of modern first principles 
pseudopotentials has allowed a vast range of problems 
to be solved for realistic materials. However, whatever 
the successes of this method, it is still unclear how non- 
local Hamiltonians should be coupled to arbitrary mag- 
netic fields. In the specific case of uniform magnetic fields 
described with the symmetric-gauge vector-potential, we 
have derived the correct Hamiltonian by developing what 
we called the gauge including projector augmented wave 
(GIPAW) method In this Letter we derive the pseu- 
dopotential Hamiltonian that describes the coupling be- 
tween the electrons and a smooth non-uniform magnetic 
field represented by an arbitrary vector-potential gauge. 

We are not alone in this quest. Recently Ismail-Beigi, 
Chang and Louie 2], to whom we refer the reader for a 
more complete summary of earlier work in this field, pro- 
posed a scheme (hereafter referred to as the ICL method) 
which sought to finally provide a rigorous derivation of 
a closed form for the coupling of nonlocal systems to ar- 
bitrary electromagnetic fields. While the widespread use 
of nonlocal pseudopotentials clearly provided the motiva- 
tion and applications for their method, they attempted 
to tackle the problem of a general nonlocal Hamiltonian. 
In contrast, this Letter focuses exclusively on the class of 
nonlocal Hamiltonians that arise due to the use of first- 
principles nonlocal pseudopotentials for electronic struc- 



ture calculations. 

In the pseudopotential approach, and in the absence 
of a magnetic field, the all-electron (AE) Hamiltonian 
H^^ = /2 + V{r), is replaced by its pseudo (PS) equiv- 
alent, TfPS = p2/2 + v\y) + Er^r • ^'(r) is a locah 
potential, and is a nonlocal operator which acts only 
within the core region of the atomic sites, at R: 



Vi = j d\'d\"\v'){v"\V^{v\v" 



(1) 



By construction, in the valence-energy-range and to 
within some controllable error, (i) the eigenvalues of H^^ 
coincide with those of H^^ , and (ii) the eigenstates of 
H^^ coincide with those of H^^ outside the core regions. 

Turning on a magnetic field, B(r) = V x A(r), the AE 
Hamiltonian becomes: 



-A(r) 

c 



V{v) 



(2) 



The question that we aim to answer in this Letter is: 
what is the PS Hamiltonian that, in presence of a mag- 
netic field, satisfies the requirements (i) and (ii) stated 
above? The PS Hamiltonian can be written as: 



Hi' 



V\r)+J2v^' + AHj^, (3) 



R 



where Ai/A remains to be determined. AH^ cannot 
be zero for all A(r), as H^^ must be gauge invariant to 
satisfy (i) and (ii). Its eigenvalues must not depend on 
the arbitrary choice of the gauge of A(r). However, as 
we will demonstrate, demanding gauge invariance alone 
is not sufficient to uniquely determine Aif^- 

To obtain H^^ we use Blochl's projector augmented 
wave (PAW) theory Q . In the PAW approach the pseud- 
isation procedure is defined as a linear transformation be- 
tween Hilbert spaces - those of the AE valence wavefunc- 
tions, and the computationally convenient PS wavefunc- 
tions. This transformation, = T|\&), can be applied 
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to obtain PS operators, O, which correspond to their AE 
counterparts, O: 



o = r^cr = o + c 



(4) 
(5) 



2J 



where we have adopted Blochl's tilde to denote a PS 
quantity. By construction, the expectation values of 
O between PS wavefunctions are equal those of O be- 
tween the corresponding AE wavefunctions. The {(jn) and 
are atomic all electron and pseudo partial waves re- 
spectively. The projector functions \pi) act only within 
some augmentation region - or core radius in the lan- 
guage of pseudopotential theory. The AE and PS par- 
tial wave coincide outside this augmentation region, and 



{P^\4'^ 



If the norms of the AE and PS partial 



waves are equal and there is just one partial wave in 
each angular momentum channel, which we shall assume 
in the following, then taking Eq. |@J with O = H^^ and 
B(r) = 0, we obtain the H^^ of norm- conserving pseu- 
dopotentials in the Kleinman-Bylander form 0|. 

In principle, in order to compute expectation values of 
physical observables from PS wavefunctions one should 
use the PS operator O. In practice, since the AE and 
PS partial waves are identical outside the augmenta- 
tion regions and have the same norm inside, the term 
C in Eq. Q can be neglected for a large class of oper- 
ators. For example, the AE perturbation Hamiltonian 
describing the coupling with a uniform electric field E is 
V^^{r) = r • E. In norm-conserving PS calculations of 
the electric-field response-properties, the C term is always 
neglected, with V^^(r) being used as the approximate 
PS perturbation potential |EMI3l- The C term can also 
be neglected for operators for which the weight is con- 
centrated away from the atoms or are almost constant 
in the augmentation regions. This is the case for oper- 
ators of the form [c • (r — R)]" in systems with a single 
PS atom centered at R. Here, and in the following, c is 
a constant vector and n > 0. In addition, the C term is 
exactly zero for the R-centered angular-momentum op- 
erator Lr = (r — R) X p, as the partial waves are eigen- 
states of |LrP and of (Lr)^. Following this reasoning 
the C term is also negligible for operators of the form 
Lr[c • (r — R)]". In contrast, C can not be neglected for 
the kinetic energy operator, or for operators that 

are general functions of the p and r operators. Indeed, by 
construction, the PS wavefunctions can be expanded on 
a much smaller set of Fourier components in the momen- 
tum space (where p is diagonal) than the corresponding 
AE wavefunctions. 

In general, because of the presence of the p operator, 
the C term corresponding to the AE magnetic perturba- 



tion Hamihonian H^^ - H^^, 



p ■ A(r) + A(r) 
2^ 



p A(r) 
2c2 



(6) 



can not be neglected, even if B(r) and A(r) are smooth 
on the scale of the core augmentation regions. However, 
as we shall show, if B(r) is smooth, and we consider a 
system with a single PS atom centered at R, there is a 
special vector potential gauge, A'(r), in which the C term 
can be neglected. This A'(r) potential can be defined 
in terms of the Fourier components be of the magnetic 
field, 



|G|<G„ 



B(r)= ^ 



(7) 



G 



where bo • G = 0, and, if the field is smooth on the scale 
of the core radius Tc, then r'cGmax ^1- In particular. 



G|<G„ 



A'(r)= E 



1 



be X (r-R)/[iG-(r-R)], (8) 



G 



with f{x) defined by, 
fix) = 2- 



l + -x + 0(x2). (9) 



Note that if the field is smooth, in the core region x ^ 1 
and f{x) can be expanded in powers of x. Since, for 
a uniform magnetic field A'(r) reduces to the symmet- 
ric gauge, A'(r) can be seen as a generalization of the 
symmetric gauge to no n- uniform magnetic fields. 

The magnetic coupling Hamiltonian in the A'(r) gauge 

is: 



AH 



AE 



|Gt<G„ 

E 

G 



bG-LR/[zG-(r-R)]+cc , A'irf 



4c 



2c2 ' 



(10) 

where cc stands for the complex conjugate. All the terms 
in the asymptotic expansion of LR/[a;] arc of the form 
Lr[c • (r - R)]" with n > 0. As a result, if B(r) is 
smooth on the scale of the augmentation region, then C 
term arising from the AH^ operator can be neglected. 
Thus the total PS Hamiltonian in the special A'(r) gauge 
is: 



1 



P + -A'(r) 

c 



1 2 



+ V\r) + V^\ (11) 



From this result we can show that the Hamiltonian for 
an arbitrary gauge A(r) is: 



H 



PS 



1 



p+ -A(r) 



+ V\r) + / d-^r'd^r"|r')(r"| x 



F^'(r',r")e"^''^— (12) 

where r r' indicates a straight line path from point 
r to point r'. To prove Eq. H12|l we have just to notice 
that is gauge invariant and reduces to Eq. Hll|l for 
A(r) = A'(r). 
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Finally, the linearity of the Hamiltonian with respect 
to the electron-ion potential can be exploited to obtain 
when there is more the one PS atom: 



Hi' = 



-A(r) 



+ V\r) + 



d\'d\" \r'){r"\ X 



(13) 



R 



We refer to as the GIPAW Hamiltonian for an arbi- 
trary magnetic field, since reduces to the Hamilto- 
nian that we have derived in our earlier work , if the 
magnetic field is uniform and if the symmetric gauge is 
used. Our new Hamiltonian holds if the magnetic field 
varies smoothly over the core region, and if the potentials 
are norm conserving. If the field varies more rapidly, or 
the norm conservation is relaxed, the C terms must be 
included in the Hamiltonian. 

Our result differs from that of ICL. The Hamiltonian 
derived in Ref. can be obtained from Eq. (|13|1 if one 
replaces our dog-leg path r' — > R — > r" with a straight 
line path r' r". The ICL Hamiltonian is also gauge 
invariant. In the presence of a magnetic field the value 
of the integral depends on the path, since V x A(r) ^ 
0. The two gauge invariant Hamiltonians are therefore 
different. 

To clarify the situation, we examine the consequences 
of these differences for the calculation of the magnetic 
susceptibility. To compare the results obtained within 
the two methods with all-electron calculations, we re- 
strict ourself to molecular systems for which it is possi- 
ble to compute the susceptibility with quantum chemi- 
cal approaches. The macroscopic magnetic susceptibility 
tensor x is defined from the second derivative of the sys- 
tem energy with respect to the external uniform magnetic 
field B: 



B • X • B = -2£;(2) 



(14) 



where i?^^^ is the second order variation of the energy 
with respect to the magnetic field: 



^ 2^[(*(0)|i/«e(eo)^('^|*(°)) 

O 



(15) 



1^^°^) and ei are the unperturbed eigenstates and eigen- 
values, g{eo) = Ee |*i°^)(^i°V(eo - Ce) and the o and 
e sums run over occupied and empty orbitals. 

In a uniform magnetic field, with the gauge A(r) ~ 
B X r/2, our GIPAW Hamiltonian, Eq. ifT^ . gives rise to 
the following perturbation Hamiltonians: 



f>(i) 

^ GIPAW 




B, 



(16) 



TABLE I: Gauge invariance test. The magnetic susceptibility 
of valence electrons (in cm'^/mole) of CH4 is calculated using 
the GIPAW and ICL approaches. Tr{x)/3 is reported as a 
function of the distance d (in a.u.) of the carbon nucleus from 
the gauge origin. The results are decomposed in terms of the 
2 contributions present in Eq. (1151 . 



d 


f/(2) hWqhW 
AGIPAW AGIPAW 


Xgipaw 


„i2) 
XlCL 


AICL 


XICL 


0.0 


-28.4 


8.6 


-19.8 


-28.4 


8.4 


-20.0 


2.5 


-68.0 


48.2 


-19.8 


-68.0 


48.0 


-20.0 


5.0 


-186.8 


167.1 


-19.8 


-186.8 


166.9 


-20.0 


7.5 


-384.9 


365.2 


-19.8 


-384.9 


364.9 


-20.0 


10.0 


-662.2 


642.5 


-19.8 


-662.2 


642.3 


-20.0 



and 



(2) 

GIPAW 



1 

8^ 



(B X r)2 + ^ B X R • ■ B X R 



R 



(17) 

'[ra, [rf3,V^^]], and a 



where = [r,V^']/i. D^^^^^ 
and /3 are Cartesian indexes. 

The corresponding perturbation Hamiltonians ob- 
tained following the ICL approach are 



(1) _ 



1 



ICL = ^ ^ ^) • ^' 



(2) 
ICL 



1 

8? 



(B X r)^ + ^ B X r • Dr • B X r 



R 



(18) 



, (19) 



where v = p -I- Er Vr. 

We compute % in molecules with both the GIPAW 
and the ICL approaches. We describe the electronic 
structure with density functional theory in the local den- 
sity approximation. We use a large-cubic-periodic su- 
percell of 6000 au'^, in order to avoid the interaction of 
the molecules with their periodic replica, and TrouUier- 
Martins pseudopotentials in the Kleinman-Bylander 
form j^. We expand the wavefunctions in plane- waves 
with a cutoff of 100 Ry. The position operator is not 
defined within periodic boundary conditions. We treat it 
approximately by constructing a periodic saw-tooth like 
function centered on the molecules. For large cells this 
operator well approximates the position operator where 
the electron density is not negligible The contribu- 
tion of core electrons to magnetic properties can not be 
neglected. This contribution is however rigid, i.e. inde- 
pendent of the chemical environment 0, . We compute 
the core contribution with an atomic code. 

Both the GIPAW and the ICL approaches are, by con- 
struction, gauge invariant. To verify that our numerical 
implementation and the use of a finite basis set preserve 
this property, we compute x for & CH4 molecule as a 
function of the distance between the gauge origin and 
the molecular center. The results summarized in Ta- 
ble|l]show that the calculated x is indeed gauge invariant. 
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TABLE II: The principal values of the magnetic susceptibil- 
ity tensor (in cm^/mole) calculated for a selection of small 
molecules. The IGAIM results are obtained from the Gaus- 
sian98 code using the aug-cc-pVQZ basis set for H and 
the aug-cc-pCVQZ basis sets for the remaining elements [TH . 
We include in the GIPAW and ICL results the rigid core con- 
tributions obtained using an atomic code. The GIPAW and 
ICL results are indicated as deviations from the IGAIM re- 
sults. The root mean square (RMS) and maximum absolute 
deviations and percentage deviations for the two methods are 
reported at the foot of the table. 

Mol. IGAIM AGIPAW AICL 



consistently closer to the all electron results than those 
calculated using the ICL method by roughly an order of 
magnitude. E.g., the ICL results for P2 deviate by almost 
one third of the total, while the maximum fractional de- 
viation for our method is less than two percent. 

We have derived, and demonstrated the practical util- 
ity of, a theory for the coupling of nonlocal pseudopo- 
tentials to arbitrary electromagnetic fields. While the 
ICL method may be the best possible for an arbitrary 
nonlocal Hamiltonian we have shown that by focusing on 
the nonlocal Hamiltonian most frequently encountered 
in electronic structure calculations we are able to make 
considerable improvements. This is of great importance 
if the results are to be compared quantitatively with ex- 
periment. 

The work of CJP was supported by the Universites 
Paris 6 and 7 in France, and the EPSRC in the United 
Kingdom. We thank the High Performance Comput- 
ing Facihty of the University of Cambridge for access to 
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IDRIS supercomputing center. We would also like to 
thank Thomas Gregor for the molecular geometries, and 
Kirk Peterson for his kind help with the Gaussian basis 
functions. 
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32.1 







while the individual terms due to H^^^ and H^"^^ clearly 
are not. 

In order to obtain an independent assessment of the 
accuracy of the two methods, we compare to a truly 
all electron method, the individual gauges for atoms in 
molecules (IGAIM) method 12j as implemented in Gaus- 
sian98 The magnetic susceptibility converges only 

slowly with Gaussian or atomic basis sets. The aug-cc- 
p(C)VxZ basis set series has been previously shown 
exhibit reliable convergence for magnetic response prop- 
erties 1^ ll^ , and we confirm this by converging the mag- 
netic susceptibility for CH4 by using up to the aug-cc- 
pCV5Z basis for C (and aug-cc-pV5Z for H). For the 
remaining calculations we use the corresponding quadru- 
ple zeta basis sets, at which level the CH4 result is con- 
verged to better than 0.1 cm'^/mole. Indeed, calculations 
using the largest basis sets rapidly become intractable for 
even moderately sized molecules. The results are sum- 
marized in Table |nl and show that the GIPAW method 
results in values for the magnetic susceptibility that are 
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